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Abstract 

We consider the force-free magnetosphere of an aligned rotator. Its 
structure is well described by a Grad-Shafranov equation for poloidal mag- 
netic field, the so-called "pulsar equation". We solve this equation by a 
multigrid method with high numerical resolution. On the fine numeri- 
, cal grid current layer along the last closed field line could be accurately 

£ — ■ incorporated into the numerical procedure and all physical properties of 

(~^) the solution such as Goldreich- Julian charge density, drift velocity, en- 

l/") , ergy losses etc. could be accurately calculated. Here we report results of 

■ the simulations. Among other interesting properties, the solution is not 

unique, i.e. the position of the last closed field line, and hence the pulsar 
f*) ■ energy losses, are not determined by the global magnetosphere structure 

I ' and depend on kinetic of electromagnetic cascades. We discuss the prop- 

erties of the solutions and their implications for pulsars. A model for 
pulsar magnetosphere evolution is proposed. In the frame of this model 
values of the pulsar braking index less than 3 have natural explanation. 



1 Pulsar Equation 

o5 ' We consider the force- free magnetosphere of an aligned rotato r. Such mag- 

netosphere is well described by the so-called p ulsar equation l|Michei Il973t 
IScharlemann fc Waeoneil Il973t lOkamotol Il974l) . For young pulsars the differ- 
ential rotation of the op en field lines could be neglected and the pulsar equation 
takes the form (see e.g. lGoodwin et aTT . 12004^1 

( X 2 - l)(d xx lP + 8 zz tp) + —- W-S^r =0 (1) 

x dip 

All coordinates are normalized to the light cylinder radius i?LC = c/£l. We 
use normalization for ip, where the dipole magnetic field function near the NS 
surface is given by 

2 

= (z 2 + Z 2 )3/ 2 ' (2) 
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The normalized poloidal current function S is given by S = {^n / c){Bj^ c / , 
where /i = BqR^ s /2 is the magnetic momentum of the NS. Magnetic field is 
expressed through the functions / and ^ = ((J./Ri,c)ip as 

B = X^i + ^I e , (3) 

vo cm 

The poloidal current S is found from the condition at the light cylinder 

a^| LC = \ ss\ c (4) 

We assume existence of an equatorial current sheet with the return current. 
The equatorial current sheet can be excluded from the numerical treatment 
by setting appropriate boundary conditions at the equatorial plane. For the 
current sheet at the boundary between closed and open field line regions this 
is not possible, and we smear the return current flowing along the last closed 
magnetic field line over the region [?/>iast — dip, V'last], see Fig. ^ 

A solution of the pulsar equatio n CI) for dipolar magnet ic field geometry 
had been found for the first time bv IContopoulos et al.1 <)l999h . They assumed, 
that t he last closed field line is fixed at the light cylinder. Later IGoodwin et alJ 
l)2004[) have found solutions for different positions xq of the point, where the last 
closed field line intersects the equatorial plane, inside the light cylinder, cc < 
1. However they have made a rather unnatural assumption, that the plasma 
pressure inside the closed field line domain is finite and plays an important role 
in the force balance across magnetic field lines. 

We perform calculations for different xq < 1 with high numerical resolution 
using grids wit h number of po i nts in each directions ranging from 2000 to 6000. 
In contrast to IGoodwin et alJ l|2004h we assume zero pressure inside the closed 
field line region, i.e. the plasma there is cold. Equation (JJJ is solved numerically 
in the domain xns < x < x ma-x , ^ns < z < z max . Boundary conditions arc shown 
in F ig. U The equation Q i s solved by a full multigrid (V-cycles) FAS scheme 
Csee iTrottenberg et all [2001). We use Gauss-Seidel smoother and SOR Gauss- 
Seidel solver at the coarsest level. Independence of the obtained results on the 
domain sizes (:ens> x max , 2ns, Zmax), dip and algorithm-specific parameters has 
been verified. 

We have expanded ip at the light cylinder in the Taylor series up to the second 
order terms, and substituted this expansion into the pulsar equation imposing 
the continuity of ip at the light cylinder. In this way we obtained an equation for 
ip at the light cylinder, accurate up to the second order term, what is sufficient 
for our numerical treatment, because we use a second order scheme. Eq.Q is 
used for determinatio n of the poloidal current. Such techniques was used by 
IGoodwin et al.l ((2004) in their simulations. It allows numerical treatment of the 
problem in a single domain, rather than doing extremely CPU time consuming 
matching of the solutions inside an outside the light cylinder, as it was made by 
IContopoulos et, alJ l(l999h . 
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2 Main Results 



• Calculations have been performed for the following values of xq- 0.15; 0.2; 
0.3; 0.4; 0.5; 0.6; 0.7; 0.8; 0.9; 0.95; 0.992. An unique solution has been 
found for each of the above Xo's. As representative examples, in Figs. 13181 
properties of solutions for three cases (xo = 0.2; 0.7; 0.99) are shown. 

• Each solution has been checked for applicability of the force-free condition 
E < B. In none of them this condition is violated, at least up to 16 light 
cylinder radii from the NS. The force-free condition can be reformulated 
as the condition on the drift velocity 



1ExB tes 

UD =c-B^ <C - (5) 



In Figs. I3I5I7I maps of u-q for three values of xo:0.2,0.6 and .99, are shown. 
One can see, that the critical value of «d = 1 is never achieved. In all 
other cases the situation is similar. 

For Xo close to 1 (xo = 0.95; 0.99; 0.992) the magnetic field strength B at 
the equatorial plane increases very rapidly when x approaches xq. The 
maximum value of B, achieved in the closed field line domain near the 
point xq, increases with increasing of xo and decreasing of dip, see F ig. U 
We in terpret this as a result of lim x _,i B = oo, predicted by IGruzinovl 
(200l). 



Energy losses of the pulsar for each Xo have been calculated numerically. 
They increase with decreasing of xo and can be described by the formula 



W = A(x ) * W md = A(x ) * HElL ( 6 ) 

o c 



where W m d are magnetodipolar energy losses. The function A(xq) can be 
fitted surprisingly well by the power law 

A(x ) = 1.41 x- 2 - 065 , (7) 

see Fig.^J We note, that analytical estimations of the energy losses based 
on the Michel poloidal current distribution for the dipolar magnetic field 
gives 

A(x ) = Xq 2 (8) 

• The total energy of the electromagnetic field in the magnetosphere E to t = 
J B ~g 7r E dV decreases with increasing of x , see Fig. ^3 

• In configurations with xo > 0.6 there is a volume return current flowing 
along the open magnetic field lines, which carries however only a small 
part of the whole return current. For xo < 0.6 the return current flows 
only along the last open field line. On Fig llll the volume current density 
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along open field lines in the polar cap of pulsar is shown. One can see, 
that the current density (in units of jcj) close to the polar cap boundary 
increases with increasing of xq. The current density does not exceed the 
Goldreich- Julian current density j GJ = p G3 c and at most field lines it is 
less than j GJ 

3 Evolution of the pulsar magnetosphere 

For each xq a special unique volume current distribution is necessary in order 
to support a force-free configuration of the magnetosphere. The current den- 
sity along open field lines in the polar cap of pulsar could adjust to the values 
required by the global magnetosphere structure (both smaller and greater than 
j G .i) when there is a particle inflow from the magnetosphere into the polar cap 
acceleration zone ijLvubarskiiL Il992) . As a source of particles in the magneto- 
sphere the outer gap cascade could be considered. On the other hand, when 
the polar cap cascade is strong and produces a lot of secondary particles with a 
wide energy distribution, some of these particles with the smallest energies, can 
be reversed in the outer magnetosphere either as a result of momentum redis- 
tribution in the turbulent outflowing plasma, or by a weak electric field arising 
as the magnetosphere will try to achieve the most energetically favorably, i.e. 
force-free, configuration. The polar cap and the outer gap cascades can not 
supply any arbitrary particle density, and, hence, can not support any arbitrary 
current density in the magnetosphere. However, a current with density greater 
or smaller than j G] can flow through the polar cap acceleration zone only if there 
is a particle inflow from the magnetosphere. The large the deviation of j from 
j G j, the larger particle flux from the magnetosphere is necessary. So, when the 
pulsar becomes older, the range of current densities, which could be supported 
by the cascades becomes smaller and the magnetosphere should change config- 
uration to a new one, where the required current density could be supported by 
the weaker cascades. 

As the electromagnetic energy of the magnetosphere decreases with increas- 
ing of xq 7 the most preferably configuration is that, where the null point (point 
where the last closed field line intersects the equatorial plane) lies at the light 
cylinder. However, in such configuration the deviation of the poloidal current 
density at the open field lines from the Goldreich-Julian current density is the 
largest and, hence, this requires a powerful source of the particles in the mag- 
netosphere. When the pulsar becomes older and cannot support such poloidal 
current, the null point moves to a new position, closer to the NS, corresponding 
to a magnetosphere configuration, where the deviation of the current density 
from j G3 is smaller and such current density could be supported by the weaker 
cascades. In the new configuration the energy losses (for the same f2) are larger 
than in the configuration with xq = 1 and, hence, the energy losses of pulsar 
decrease with time slower than it is given by the magnetodipolar formula. If at 
some time the dependence of xq on the pulsar age and, hence, on Q is given by 
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x oc Qt, then using eqs. 0,® we get 

W(xn a , a = 4- 2.065 £, (9) 

£ is a complicated function of the pulsar age i, but as xq decreases with t (see 
above), it is always positive, £(t) > 0. So, the braking index 

fin , x 

n= — =a-l = 3 - 2.065 £, (10) 

is always less than 3!. 

Detailed description of the numerical method and results of calc ulation, as 
well a s discussion of their implication for pulsar physics, are given in iTimokhirJ 
l|2005h . 

This work was supported by RFBR grant 04-02-16720. 
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Figure 1: Calculation domain and boundary conditions 
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Figure 2: Magnetic field strength at the equatorial plane as a function of x for 
different xq and widths of the current sheet. 
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Figure 3: Magnetic field lines configuration (radial black lines) for xq = 
0.99, dip = 0.08. The last closed field line is shown by the thick red line. Con- 
tours of the absolute values of the drift velocity normalized to the speed of light, 
|ud|/c, are shown by the blue vertical lines. 




0.5 1 1.5 



Figure 4: Central parts of the magnetosphcre for xq = 0.992, dip = 0.04. Nota- 
tions similar to ones in Fig. Also line where p G3 changes sign and j po i = 
are shown. 
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Figure 5: Drift velocity and magnetic field configuration for xq — 0.7, dip = 0.03 




Figure 6: Inner parts of Fig. 
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Figure 7: Drift velocity and magnetic field configuration for xq = 0.2, dip = 0.03 





Figure 10: Total electromagnetic energy of two different volumes as a function 
of xq . The values were normalized to the energy in the corresponding volumes 
for Xq — 0.15. 
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Figure 11: Poloidal current density in the polar cap of pulsar for different val- 
ues of xo, normalized to the corresponding Goldreich- Julian current density 
as a function of the colatitude, normalized to the colatitude of the polar cap 
boundary 8 pc . Michel current density is shown by the dashed line 
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